The Fish Pathogen “Candidatus Clavichlamydia salmonicola”—A Missing Link in the Evolution of Chlamydial Pathogens of Humans

Abstract Chlamydiae like Chlamydia trachomatis and Chlamydia psittaci are well-known human and animal pathogens. Yet, the chlamydiae are a much larger group of evolutionary ancient obligate intracellular bacteria that includes predominantly symbionts of protists and diverse animals. This makes them ideal model organisms to study evolutionary transitions from symbionts in microbial eukaryotes to pathogens of humans. To this end, comparative genome analysis has served as an important tool. Genome sequence data for many chlamydial lineages are, however, still lacking, hampering our understanding of their evolutionary history. Here, we determined the first high-quality draft genome sequence of the fish pathogen “Candidatus Clavichlamydia salmonicola”, representing a separate genus within the human and animal pathogenic Chlamydiaceae. The “Ca. Clavichlamydia salmonicola” genome harbors genes that so far have been exclusively found in Chlamydia species suggesting that basic mechanisms important for the interaction with chordate hosts have evolved stepwise in the history of chlamydiae. Thus, the genome sequence of “Ca. Clavichlamydia salmonicola” allows to constrain candidate genes to further understand the evolution of chlamydial virulence mechanisms required to infect mammals.


Introduction
Chlamydiae are a large phylum of obligate intracellular bacteria infecting a broad spectrum of eukaryotic hosts ranging from protists to invertebrates and mammals (Taylor-Brown et al. 2015;Collingro et al. 2020). Members of six chlamydial families, namely the Piscichlamydiaceae, Parilichlamydiaceae, Clavichlamydiaceae/Chlamydiaceae, Parachlamydiaceae, Simkaniaceae, and Rhabdochlamydiaceae are able to infect a variety of fish, where they are-together with a few other bacteria-associated with the poorly understood gill and skin disease epitheliocystis (Toenshoff et al. 2012;Stride et al. 2014;Katharios et al. 2015;Pawlikowska-Warych and Deptuła 2016;Seth-Smith et al. 2016;Blandford et al. 2018). Members of the Piscichlamydiaceae, Parilichlamydiaceae, and Clavichlamydia are particularly prevalent. They seem to be restricted to fish as they have not been detected in other hosts so far. Any attempts to culture these bacteria in surrogate host systems in the lab have failed so far (Blandford et al. 2018). "Candidatus Clavichlamydia salmonicola" (hereafter referred to as Cl. salmonicola), infects farmed and wild salmonids in freshwater and often co-occurs with Piscichlamydia species (Rourke et al. 1984;Bradley et al. 1988;Karlsen et al. 2008;Mitchell et al. 2010;Schmidt-Posthaus et al. 2012;Guevara Soto et al. 2016a, 2016b, 2017.
Candidatus Cl. salmonicola represents an early branching genus within the Chlamydiaceae, which includes important human and animal pathogens like Chlamydia trachomatis, Chlamydia pneumoniae, and Chlamydia psittaci (Elwell et al. 2016). Since attempts to cultivate Cl. salmonicola in fish and insect cell lines as well as amebae were not successful, and as publicly available metagenomic studies have not captured clavichlamydiae so far (Köstlbacher et al. 2021b), we here applied a cultivation-independent, targeted metagenomic approach to obtain insights into the genetic repertoire of this chlamydial species. In order to study common chlamydial virulence mechanisms in vertebrate hosts, that is to identify genes involved in host adaptation, host-specificity, and virulence, we then compared the Cl. salmonicola genome with other chlamydial genomes with a focus on the closely related members of the Chlamydiaceae and the recently proposed family Sororchlamydiaceae (Dharamshi et al. 2020(Dharamshi et al. , 2022 as well as the only other fish pathogens with available genome sequences, members of the Parilichlamydiaceae (Taylor-Brown et al. 2017, 2018.

Results and Discussion
General Genome Features Total genomic DNA was extracted from infected gill tissue of Atlantic salmon (Salmo salar) from a commercial fish farm in the northwest of Ireland in February 2009 (Mitchell et al. 2010) and sequenced using a combination of pyrosequencing and Illumina technologies.
The assembly of the Cl. salmonicola genome consists of 28 scaffolds with a total length of 1,392,990 bp (25 × coverage) and a G+C content of 32.48% (supplementary table S1, Supplementary Material online). According to analysis with CheckM2 v0.1.3 (Chklovski et al. 2022), the genome assembly has a completeness of 94.68% with 0.27% contamination, suggesting a theoretical genome size for Cl. salmonicola of roughly 1.47 Mbp. The assembly includes a complete rRNA gene operon on the largest 470,378 nt scaffold, as well as 37 tRNA genes, and 1,105 protein-coding sequences (CDS). One of the scaffolds (contig 16) represents a plasmid of 8,040 bp (9 CDS, G+C content 28%), which is highly similar to the small plasmids of other Chlamydiaceae species (supplementary table S1, Supplementary Material online, supplementary fig. S1, Supplementary Material online) (Köstlbacher et al. 2021a).

Clavichlamydia Salmonicola Represents a Novel Genus Within the Chlamydiaceae
Previous 16S rRNA gene-based phylogenies consistently suggested Cl. salmonicola as representative of the familylevel lineage Clavichlamydiaceae and being the closest relative of the human and animal pathogenic Chlamydiaceae (Karlsen et al. 2008;Stride et al. 2014 (Vorimore et al. 2021). While ANI serve well for species delineation, AAI are better suited to resolve genus-and family-level clades (Konstantinidis et al. 2017 (Dharamshi et al. 2020;Köstlbacher et al. 2021b). In a new study, four additional sponge-associated MAGs were published, and the name Sororchlamydiaceae was proposed for this clade (Dharamshi et al. 2022). Phylogenetic analyses in our study and comprehensive phylogenetic analyses of the chlamydial phylum in a recently published study further corroborate the phylogenetic placement of Clavichlamydia close to Chlamydia (Dharamshi et al. 2022(Dharamshi et al. , 2023. Phylogenies inferred for the 16S rRNA gene as well as the single copy marker genes show that sponge-associated chlamydiae are basal to fish (Clavichlamydia), amphibian (Amphibiichlamydia, only represented by 16S rRNA gene data so far), bird (Chlamydiifrater), and amphibia/bird/reptile/mammal (Chlamydia) infecting chlamydial genera (

Virulence-Associated Genes Shared Amongst Chlamydial Lineages
The obligate intracellular chlamydial lifestyle is ancient and shared by all known chlamydiae (Collingro et al. 2020). Their biphasic developmental cycle alternates between an extracellular infectious stage, the elementary body, and an intracellular replicative form, the reticulate body (Elwell et al. 2016). Chlamydiae employ a type III secretion (T3S) system, which has already been present in the last common chlamydial ancestor, and effector proteins as the main toolbox to infect and interact with their eukaryotic host cells (Peters et al. 2007;Mueller et al. 2014;Dharamshi et al. 2023). After host cell entry, chlamydiae usually reside within a host-derived vacuole called inclusion (Elwell et al. 2016). Various chlamydial proteins are present in the inclusion lumen or inserted into its membrane (Bugalhão and Mota 2019; Gitsels et al. 2019).
Our analysis of orthologous groups (OGs) of proteins in 49 chlamydial genomes confirmed the presence of wellconserved genes indicative of the basic host-dependent chlamydial lifestyle in Cl. salmonicola ( fig. 2) (Köstlbacher et al. 2021b). These genetic traits include the T3S system and various effector proteins, genes for host cell adhesion, acquisition of host nutrients, proteases and kinases to modulate the host cell, and major transcriptional regulators of the developmental cycle ( fig. 2).
Chlamydia species encode a number of additional gene families that so far have not been detected in any other chlamydiae and have been associated with virulence in humans and animals. With the availability of MAGs from Sororchlamydiaceae and genomes of the novel genera Clavichlamydia and Chlamydiifrater it is now possible to analyze differential presence of these genes within the Chlamydiaceae and their sister family. Interestingly, 23 of these gene families so far exclusively present in Chlamydiaceae are also found in Sororchlamydiaceae genomes ( fig. 2). Among them are some genes that have previously been shown in Chlamydia species to be important for host cell interaction. This includes the T3S system needle length determinator CdsP (Lorenzini et al. 2010), and the chaperones CopB2, CopD2, and CopD1, which are located at the tip of the T3S system needle and considered to be essential for infection of nonphagocytic host cells (Matteï et al. 2011;Bulir et al. 2014). Some other T3S effectors (CT_195 and CT_656) some of which interact with the host cell's endosomal sorting complexes required for transport (ESCRT) machinery (CT_619, CT_712; Vromman et al. 2016) as well as inclusion proteins involved in the acquisition of host cell lipids (Cap1) are shared too ( fig. 2) (Bugalhão and Mota 2019). The transcription of many genes expressed late during the developmental cycle is regulated by σ 28 in Chlamydia species (Yu et al. 2006). This transcription factor is also encoded in spongeassociated Sororchlamydiaceae and all Chlamydiaceae genomes, whereas it is absent in other chlamydiae ( fig. 2) (Domman and Horn 2015). A recent study modeling ancestral chlamydial genomes and their evolution showed that the above-mentioned genes with the exception of σ 28 appeared first in the last common ancestor of Sororchlamydiaceae and Chlamydiaceae, but their origin is unknown (Dharamshi et al. 2023) (supplementary table S4, Supplementary Material online). The acquisition of these genes might have primed members of these chlamydial lineages for the infection of metazoan hosts.
Gene families exclusively present in all Chlamydiaceae (n = 27) include additional T3S genes like copB1 and T3S effectors important for the entry of nonphagocytic host cells (tarP) (Jewett et al. 2006;Ghosh et al. 2020), potentially modulating the host cell cycle (CT_847; Elwell et al. 2016), host lipid acquisition (CT_618), or inclusion proteins (CT_006, CT_058, CT_440, incC; Bugalhão and Mota 2019) ( fig. 2). Again, apart from tarP (see below), these gene families have first occurred in the ancestor of Chlamydiaceae, and their initial source is unknown as there are no known homologs in current sequence databases (Dharamshi et al. 2023 Chlamydiaceae also harbor a small plasmid (8-9 kb), which is highly conserved in gene synteny (Supplementary fig. 1, Supplementary Material online). The Cl. salmonicola plasmid encodes all eight genes known from Chlamydia plasmids (30.1-58.9%, mean 39.9% amino acid identity) and an additional partial DNA helicase. Two of the genes on the chlamydial plasmid (pgp3, pgp4) were previously considered to be specific to members of the genus Chlamydia but are indeed found on all Chlamydiaceae plasmids. These genes have been implicated in chlamydial niche differentiation towards higher animals, and they have recently been proposed to be involved in the formation of putative outer membrane vesicles delivering proteins into the host cytosol (Köstlbacher et al. 2021a;Lei et al. 2021).
Within Chlamydiaceae, there are 57 orthologous genes exclusively shared by members of the genera Chlamydiifrater and Chlamydia ( fig. 2). This includes the membrane proteins OmcA and CT_814.1, the histone-like protein HctB, and T3S effectors affecting inclusion growth (taiP; Hamaoui et al. 2020) and interacting with the host ESCRT machinery (CT_711; Bugalhão and Mota 2019) ( fig. 2, supplementary table S4, Supplementary Material online). Again, some inclusion proteins were gained for potential host interaction and modulation including incB, incS, CT_082 and CT_083, and mrcA promoting inclusion extrusion and interfering with host cell homeostasis (Nguyen et al. 2018;Chamberlain et al. 2022). Of note, no gene is shared solely between Chlamydiifrater and Clavichlamydia, and only two genes are exclusively present in Chlamydia and Clavichlamydia genomes. One of them representing incV important for directing the endoplasmic reticulum towards chlamydial inclusions (Ende et al. 2022). Finally, of the many orthologous genes previously only known to exist in Chlamydia species still 48 are specific for members of this genus, potentially reflecting genes necessary for the interaction with bird, reptile, and mammalian cell types or tissues. Taken together, the in depth-analysis of orthologous genes shared between different chlamydial lineages, especially those shared between and amongst Sororchlamydiaceae and genera in the Chlamydiaceae, revealed major differences in genes involved in host interaction and manipulation. Our findings suggest a sequential acquisition of these genes from Sorochlamydia-like ancestors infecting basal animal clades such as sponges to Clavichlamydiae infecting fish; and further from Chlamydiifrater infecting birds to Chlamydia species infecting amphibians, reptiles, birds, and mammals.

Evolution of the Chlamydial T3S Effector translocated actin-recruiting effector
The translocated actin-recruiting effector (TarP) is a welldescribed key virulence factor necessary for the entry of nonphagocytic host cells (reviewed in Caven and Carabeo 2020), which is specific to members of the Chlamydiaceae ( fig. 2). Although TarP has no orthologs in the eggNOG database (Dharamshi et al. 2023), it forms an orthogroup in OMA orthologs (Altenhoff et al. 2021) with 28 distantly related genes from nonchlamydial organisms (supplementary table S4, Supplementary Material online). Phylogenetic analysis of this orthogroup indicates that TarP was acquired by the last common ancestor of the Chlamydiaceae ( fig. 3A). Within the family, the Clavichlamydia TarP represents the deepest branching lineage. The genera Chlamydiifrater and Chlamydia are not well-separated due to the lack of resolution in this tree, but within both genera, the gene tree is congruent with the species tree, together suggesting maintenance and further diversification of TarP within the Chlamydiaceae. It has been shown previously that TarP phylogeny in C. trachomatis correlates with tissue tropism (Lutter et al. 2010).
The recent progress of methods in genetic manipulation of chlamydiae has enabled a better understanding of TarP  and its interaction with diverse host factors (Caven and Carabeo 2020). Although the functional domains of TarP have been characterized mainly in C. trachomatis, their existence in C. pneumoniae has recently also been verified experimentally (Braun et al. 2019). Notably, these domains are absent in the closest related nonchlamydial orthologs of TarP (i.e., in Alteromonas macleodii).
A more detailed analysis of the domain architecture of chlamydial TarP revealed an interesting pattern of domain acquisition ( fig. 3B). All chlamydial TarP proteins have a proline-rich domain (PRD) and at least two G-or F-actin binding domains (ABD or FAB, respectively). PRD is responsible for TarP oligomerization in the host cell, followed by the formation of a G-actin nucleus via the concentrated localization of ABD and subsequent actin polymerization. This eventually results in the remodeling of the host cell cytoskeleton and a phagocytosis-like internalization of chlamydial elementary bodies (Caven and Carabeo 2020). While a single ABD is sufficient for host cell entry (Jewett et al. 2010;Parrett et al. 2016), members of the genus Chlamydia have adopted additional ways to recruit host actin with TarP. Their TarP proteins harbor vinculin-binding domains (VBD, n = 1-3), which account as a second path to indirectly polymerize F-actin at the  chlamydial entry site (Thwaites et al. 2015). A further step in Chlamydia adaptation to remodel actin-but missing in C. pneumoniae-is via domains rich in leucine and aspartic acid (LD-domains) located in FABs (Thwaites et al. 2014). Finally, TarP of C. trachomatis possesses a tyrosine-rich repeat region, which not only plays a role in actin polymerization but also potentially affects the host cell cycle (Shehat et al. 2021). Recently, a postinvasion function of TarP has been discovered, in which LD and VBD domains inhibit cell motility by binding focal adhesion kinase and vinculin, which leads to less cell shedding in epithelial tissue (Pedrosa et al. 2020). This could represent a mechanism of chlamydiae to maintain infections in high-turnover tissues (Pedrosa et al. 2020). Taken together, the acquisition of novel domains in TarP might reflect stepwise adaptations within Chlamydiaceae to refine their potential to manipulate their host cells.
A similar pattern of domain architecture evolution can be observed in the inclusion membrane protein IncV (shared between Clavichlamydia and Chlamydia species; fig. 2, supplementary table S4, Supplementary Material online). This protein establishes inclusion membrane contact sites (MCS) with the endoplasmic reticulum for the acquisition of host lipids (Ende et al. 2022). While the N-terminal part of the protein is well-conserved in its amino acid sequence and predicted transmembrane helices in Cl. salmonicola, the protein is shorter than in Chlamydia species and lacks all characterized motifs interacting with mammalian cells to establish the MCS-like phosphorylation sites, FFAT motifs, and VAMP-associated proteins recruitment sites (Ende et al. 2022). Other Chlamydiaceae virulence genes, however, do not differ in motifs or domains and are therefore likely functionally conserved since their origin. CT_712 and CT_619, for instance, T3S effectors interacting with the host ESCRT machinery contain an N-terminal coiled-coil domain for binding host TSG101 and a C-terminal DUF582 domain (Pfam PF04518) in all orthologs (Vromman et al. 2016). In summary, the detailed analysis of domain architecture provides valuable insights into the evolution of chlamydial effector proteins and suggests domain acquisition as a potential mechanism of sequential evolution of host interaction, and eventually host-specificity and tissue tropism.

The Chlamydial Plasticity Zone
Genomes of Chlamydia species are generally highly similar and contain only a few regions with a higher degree of variation. One of these variable regions is the plasticity zone (PZ), which ranges in different Chlamydia species from ∼5,500 to ∼55,500 bp including between six and 49 genes (Hölzer et al. 2020). The PZ is located between the genes encoding subunits of the acetyl-CoA carboxylase (accBC) and the purine biosynthesis genes guaAB, but the latter genes are absent in some chlamydial PZs (Nunes and Gomes 2014;Hölzer et al. 2020). The gene content within the PZ differs drastically between species but can contain important genes implicated in virulence and possibly tissue tropism, including a cytotoxin, phospholipase D, membrane attack complex component (MAC)/perforin, and a tryptophan biosynthesis gene cluster (Nunes and Gomes 2014) (fig. 4).
So far, the PZ has only been described in genomes of Chlamydia species. However, with the availability of genomes of members of two novel genera within Chlamydiaceae and Sororchlamydiaceae, we were able to check for its presence also in these genomes. While the accBC genes are present in members of Sororchlamydiaceae, we could not detect any other genes specific for the Chlamydia PZs and only a partial conservation of gene synteny over a large 200 kb genomic region (contig 1 nucleotides 110800-306700; fig. 4). In the clavichlamydial genome none of the cytotoxin, MAC/perforin, or tryptophan synthesis genes were detected, but there is some synteny with other (housekeeping) genes in Chlamydiifrater and the C. trachomatis PZ including genes right up-and downstream of the virulence-associated PZ genes. In comparison to Sororchlamydiaceae, the respective genomic region is condensed ( fig. 4). A copy of the phospholipase D gene outside the PZ, which is found in all chlamydial genomes, is encoded elsewhere in the clavichlamydial genome and seems to have been duplicated (CLAVI_RS02740, CLAVI_RS02745, 34.9% aa identity) (Dimond and Hefty 2021). In Chlamydiifrater, the cytotoxin CT_166 is present in the respective region, which also shows more synteny including some rearrangements to the C. trachomatis PZ ( fig. 4). Together, this suggests that a genomic region resembling the extant PZ was present in the last common ancestor of Chlamydiaceae and has diverged in members of the different genera. It is more similar in members of the sister genera Chlamydiifrater and Chlamydia, suggesting a stepwise evolution of the chlamydial PZ in the course of the adaptation to life inside fish, bird, reptile, and mammalian hosts and tissues.

Genomic Adaptations to Fish Hosts
Candidatus Cl. salmonicola encodes 175 singleton genes and 25 larger gene families (with up to 54 members) that have no or only weak homologies to any known genes. Some genes in these large clavichlamydial gene families harbor motifs similar to those present in transport-associated genes, and in 14 of these OGs, 33-100% of the genes contain predicted transmembrane helices, together suggesting functions in transport, adhesion, or cell integrity. The remaining eight larger clavichlamydial families show no homology to any known motifs or protein domains. Examining the function of genes belonging to Fish Pathogen "Candidatus Clavichlamydia salmonicola" GBE Genome Biol. Evol. 15(8) https://doi.org/10.1093/gbe/evad147 Advance Access publication 24 August 2023 these conspicuous gene families might help to better understand how Cl. salmonicola cell biology and host interaction differs from other chlamydiae. The general pattern of expanded gene families, though, has been observed previously in other chlamydial species and might increase the potential for host interaction and host range (Domman et al. 2014).
Since the only known fish-pathogenic chlamydiae with available genome sequences are members of the Parilichlamydiaceae and Cl. salmonicola, we were interested in genes shared exclusively by members of these two families. To our surprise, we were not able to detect any genes exclusively shared between genomes of Parilichlamydiaceae and Cl. salmonicola. Since all genomes available for members of those two families are draft genomes, shared genes might either have been missed in our analysis, or the genetic basis of the interaction with fish hosts is indeed different. The genomes of Parilichlamydiaceae are the smallest chlamydial genomes known to date and underwent massive loss of metabolic potential ). Interestingly, we found some congruence in genes and pathways absent in Parilichlamydiaceae and the Cl. salmonicola genomes. Both groups lack a number of genes for cell division, lipopolysaccharide, and glycerolipid biosynthesis, suggesting that mechanisms for cell division and cell wall generation might differ in these fish pathogens from other chlamydiae ( fig. 2) (Pillonel et al. 2018;Taylor-Brown et al. 2018). In addition, genes for heme, riboflavin, methylene tetrahydrofolate biosynthesis, the shikimate pathway, and the tricarboxylic acid cycle are missing in genomes from both groups . As these pathways include a large number of genes, often distributed throughout the genome, it seems unlikely that their lack is purely due to their absence in the assemblies of the (high-quality) draft genomes of Parilichlamydiaceae and Cl. salmonicola. While complete genome sequences for these species would be needed to verify these observations, current evidence suggests a common pattern regarding the absence of certain metabolic traits between chlamydiae infecting fish.
As the phylogenetic placement of Parilichlamydiaceae is not fully resolved yet (Dharamshi et al. 2023), it remains to be seen whether the potential loss of these genes is due to convergent evolution in a similar vertebrate host as suggested earlier , or if these organisms share a closer relationship than previously assumed. In both scenarios, the lack of important metabolic traits might be the reason for the restriction of these chlamydiae to fish hosts. To gain more insights into these 10 kb  important questions, successful isolation and cultivation of these chlamydiae would be of great importance.
The genome sequence of Cl. salmonicola represents the first available genome sequence of the most deeply branching genus within the Chlamydiaceae. By comparative analysis including other recently published close relatives (Vorimore et al. 2021;Dharamshi et al. 2022), we were able to identify a set of genes that seem to play a role in the adaptation of chlamydiae from less complex metazoan hosts to mammals. Starting from sponge-associated Sorochlamydiaceae to fish and bird infecting Cl. salmonicola and Chlamydiifrater species, to bird, reptile, and mammalian Chlamydia pathogens, we observed evidence for stepwise adaptations in the organization of the PZ, virulence gene content and presence of protein domains of key virulence-associated genes. The genome sequence of Cl. salmonicola thus constitutes an important puzzle piece in our understanding of chlamydial biology and evolution.

Sample Preparation
Twenty Atlantic salmons (Salmo salar) were randomly sampled by veterinarians from a commercial fish farm in the northwest of Ireland in February 2009 (Mitchell et al. 2010). Fish were held in fresh water tanks supplied with river water. All fish were from the same cohort and had a body weight of 40-60 g. The fish showed no clinical abnormalities or pathological symptoms.

DNA Extraction and Quality Control
Genomic DNA was isolated from five gill arches of eight individuals respectively. Gills were homogenized in a Dounce tissue grinder (Wheaton) in buffer A (35 mM Tris-HCl, 25 mM KCl, 10 mM MgCl 2 , 250 mM sucrose (all Sigma Aldrich), pH 7.5) containing 250 mM EDTA (Carl Roth). The suspension was filtered through a 5 µm syringe filter and centrifuged for 10 min at 6,000 rpm and 4 °C. The pellet was washed twice in buffer A, resuspended in buffer A containing 10 U DNase I (Thermofisher Scientific), and incubated for 1 h at 37 °C. DNase I digestion was stopped with 50 mM EDTA and centrifuged as above. The resulting pellet was washed in buffer A with 250 mM EDTA and resuspended in TE buffer (Thermofisher Scientific). Bacterial DNA was purified using a sodium dodecyl sulphate-based method including 1% hexadecylmethylammonium bromide (CTAB, Sigma Aldrich) and 200 µg/ml Proteinase K (Sigma Aldrich) in the extraction buffer (Zhou et al. 1996). Genomic DNA was stored at −20 °C until further use.
Semiquantitative PCR assays were performed to estimate the relative abundance of rRNA genes in the extracted genomic DNA. All PCR reactions were performed with 1 µl template DNA, 1 unit of Taq DNA polymerase, 10 x Taq buffer with KCl and 2 mM MgCl 2 , 0.2 mM of each deoxynucleotide (all Thermofisher Scientific) and 50 pmol of each primer in a total volume of 50 µl. Negative (without template) and positive controls in defined declining template concentrations were included in all semiquantitative PCR assays. The presence and size of amplicons were checked by gel electrophoresis and ethidium bromide or SYBR Green I (Sigma Aldrich) staining. Only low amounts (<10 pg/µl) of eukaryotic DNA were detected by the general 18S rRNA gene primer pair 18SF/18SR (supplementary

Genome Sequencing and Comparative Genome Analysis
Library preparation and sequencing of genomic DNA were performed at Agowa GmbH (Berlin, Germany) using a 454 GS-FLX Titanium pyrosequencing platform and at the Vienna BioCenter Core Facilities Next-Generation Sequencing Unit (https://www.viennabiocenter.org/ facilities/) using an Illumina HiSeq2000 instrument to generate paired-end reads of ∼125 bases according to standard procedures. After quality filtering with BBDuk (minlen = 50 qtrim = rl trimq = 25 ktrim = r k = 25 mink = 11) from the BBTools package (v34.24, Bushnell), the reads were used for a hybrid assembly with SPAdes v3.6.0 (Bankevich et al. 2012). The assembly was screened for completeness and contamination with CheckM2 v0.1.3 (Chklovski et al. 2022). Genome annotation was performed with ConsPred 1.22 (Weinmaier et al. 2016 All proteins encoded in these genomes were clustered into OGs with OrthoFinder 2.5.4 with default parameters (Emms and Kelly 2019). For the final analysis, 49 chlamydial genomes representing well-known chlamydial families were included (supplementary table S2, Supplementary Material online). A gene was considered as being present in a chlamydial family-level lineage, when more than 50% of the genomes of this family included in the analysis encoded the gene. The presence of genes in chlamydial genomes was visualized with the R package UpSetR 1.4.0 (Conway et al. 2017). Gene synteny in the chlamydial PZ and plasmids were visualized with the R package genoPlotR 0.8.11 (Guy et al. 2010).
The analysis of protein domains was performed by combining the comparison of ClustalW sequence alignments and InterProScan 5.62-94.0 screens at EMBL-EBI (Madeira et al. 2022).

Phylogenomic and Phylogenetic Analysis
The protein sequences of 43 conserved checkM single copy marker proteins were extracted and aligned in CheckM v1.2.2 with the "tree" workflow (supplementary table S3, Supplementary Material online) (Parks et al. 2015). We performed model testing and maximum likelihood phylogenies with IQ-TREE 2.2.2.3 (Nguyen et al. 2015) under the empirical LG model (Le and Gascuel 2008). The optimal model was determined with "-m TESTNEW" (Kalyaanamoorthy et al. 2017), including the empirical mixture models C10-C60 with the "-madd" option (best model: LG+C60+G +F) (Quang et al. 2008). We inferred 1,000 ultrafast bootstrap replicates (Hoang et al. 2017) with the "-bnni" option for bootstrap tree optimization and 1,000 replicates of the SH-like approximate likelihood ratio test (Guindon et al. 2010). The initial species tree was then used as a guide tree for posterior mean site frequency (PMSF) modeling (Wang et al. 2018) under the LG+C60+G+F model, and 100 nonparametric bootstraps were inferred. Fourteen Planctomycetota/Verrucomicrobiota species served as the outgroup for rooting the phylogenetic trees (supplementary table S2, Supplementary Material online).
As no other genome sequences are currently available for Cl. salmonicola, we performed additional phylogenetic tree inference using 16S rRNA gene sequences to analyze the affiliation of the genome sequence derived in this study with previously published Cl. salmonicola sequences (Karlsen et al. 2008;Guevara Soto et al. 2016a). We used the dataset described in Köstlbacher et al. (2021b) and added three additional clavichlamydial, two Chlamydiifrater, and one Sororchlamydiaceae 16S rRNA gene sequences (Karlsen et al. 2008;Guevara Soto et al. 2016a;Vorimore et al. 2021;Dharamshi et al. 2022). The sequences were aligned with SINA (Pruesse et al. 2012) and trimmed with trimAl "-gappyout" (Capella-Gutiérrez et al. 2009). After performing model finding with IQ-TREE 2.2.3 (Kalyaanamoorthy et al. 2017), the phylogenetic tree was inferred with the best-fitting model SYM+R10 and 100 nonparametric bootstraps.
In order to check for the presence of homologs of chlamydial virulence genes in other organisms, the dataset published in Dharamshi et al. (2023) was used. In addition, genes of interest were checked for homologs against OmaGroups in the Omabrowser (Altenhoff et al. 2021) (supplementary table S4, Supplementary Material online).
All phylogenetic trees generated in this study were visualized with iTOL version 6.7.5 (Letunic and Bork 2021).

Supplementary Material
Supplementary data are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/). number WTCQ00000000. The NCBI BioProject and BioSample numbers are PRJNA492195 and SAMN10090347, respectively.